% This script sequentially reads *.res files that are generically named after its period number
% and spherical harmonics numbers. It then interpolate the magnetic or
% electric fields at the observator/station locations specified and creates
% for each observatory and for each spherical harminc coefficient, the real
% and imaginary parts of the unit responses as a function of frequency

clear;

basedirectory = '/Users/manojnair/projects/GIC/Modelling/'

% There are two responses available. The Files with 3e3 are modeled with
% Lithospheric resitivity 3000 Ohm-m and files 3e2 are modelled with
% lithopsheric reisitivity 300 Ohm-m. The difference is maily along
% conducitivity contrast regions. It is important that the file names
% (input and output) to be consistant. 

Res_directory = 'Results_KO06_3e3';
Generic_name  = 'KO_06_gm00_lith_3e3_P'; % Example file name KO_06_gm00_lith_3e3_P14_q11.res
Output_directory_name = 'KO_06_3e3_elec';

% NHarmonics is the number of modeled harmonics. In this case it is degree
% 3 and (N*(N+2) is 15
NHarmonics = 15;

% N periods is the number of periods modeled. Here there are 50
NPeriods   = 50;

% The original model was done on 1x1 degree
nt = 180;
np = 360;

dp = 360/np; dt = 180/nt;

theta = (dt/2 : dt : 180)';
phi   = (dp/2 : dp : 360)';


[X,Y] = meshgrid(phi,theta);

harm( 1,:) = 'q10';
harm( 2,:) = 'q11';
harm( 3,:) = 's11';
harm( 4,:) = 'q20';
harm( 5,:) = 'q21';
harm( 6,:) = 's21';
harm( 7,:) = 'q22';
harm( 8,:) = 's22';
harm( 9,:) = 'q30';
harm(10,:) = 'q31';
harm(11,:) = 's31';
harm(12,:) = 'q32';
harm(13,:) = 's32';
harm(14,:) = 'q33';
harm(15,:) = 's33';
harm(16,:) = 'q40';
harm(17,:) = 'q41';
harm(18,:) = 's41';
harm(19,:) = 'q42';
harm(20,:) = 's42';
harm(21,:) = 'q43';
harm(22,:) = 's43';
harm(23,:) = 'q44';
harm(24,:) = 's44';


% RUN THIS FIRST WITH MAGNETIC DATA. DO NOT LOAD X,Y,Z data from these
% files
%load /Users/manojnair/projects/GIC/obs_apr_2000.mat fday phi_gg theta_gg name
%load /Users/manojnair/projects/GIC/obs_jul_2012.mat fday phi_gg theta_gg name

%---- new modeling -----

%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Apr_2000.mat fday phi_gg theta_gg name
%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Jul_2000.mat fday phi_gg theta_gg name
%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Oct_2003.mat fday phi_gg theta_gg name
%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Aug_2000.mat fday phi_gg theta_gg name
%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Mar_2001.mat fday phi_gg theta_gg name
%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Nov_2003.mat fday phi_gg theta_gg name
%load /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_NovSq_2003.mat fday phi_gg theta_gg name

% FOR ELECTRIC FIELD CALCULATIONS RUN THIS SECOND (AFTER THE ABOVE)
%load /Users/manojnair/projects/GIC/OBEM/obem_elec_apr00.mat fday phi_gg theta_gg names
%load /Users/manojnair/projects/GIC/obs_oct_2011.mat fday phi_gg theta_gg names
%load /Users/manojnair/projects/GIC/usarray/201110/usarray_elec_oct11.mat fday phi_gg theta_gg names
%load /Users/manojnair/data/usarray/201207/usarray_elec_jul2012.mat fday phi_gg theta_gg names
%load /Users/manojnair/projects/GIC/obs_nov_2003.mat fday phi_gg theta_gg name; % theta_gg already in co_latitudes
%load 'obsdata_elec.mat' fday phi_gg theta_gg name; % theta_gg already in co_latitudes

% ---- new modeling ------
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Apr00.mat  fday phi_gg theta_gg names
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Jul00.mat  fday phi_gg theta_gg names
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Oct03.mat   fday phi_gg theta_gg names
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Aug00.mat  fday phi_gg theta_gg names
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Mar01.mat fday phi_gg theta_gg names
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Nov03.mat fday phi_gg theta_gg names
load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_SqNov03.mat fday phi_gg theta_gg names



%convert latitude to co-latitude. NOTE. This is a place we can mess up the
%scripts. The interpolation below assumes that the latitudes are in
%co-latitudes. Do the step below only if the latitudes are not in
%co-latitudes.
%theta_gg = 90 - theta_gg;

%Rotation of the coordinate system GG->GM. It returns the geomagnetic
%latitude in co-latitude
[theta_gm, phi_gm]  = gg2gm_new(theta_gg, phi_gg);


iperiod = 1;

for iP = 14 : NPeriods
    
    for iH = 1 : NHarmonics
        
        FileName = [basedirectory Res_directory '/' Generic_name ...
            sprintf('%02d_',iP) harm(iH,:) '.res'];
        
        [Hp, Ht, Hr, Et, Ep, ...
            Period, np, nt] = read_x3dg_res(FileName);
        
        Obs_interp(1,iperiod, :, iH) = interp2(X, Y, Ht, phi_gm, theta_gm,  'linear');  %Ht
        Obs_interp(2,iperiod, :, iH) = interp2(X, Y, Hp, phi_gm, theta_gm,  'linear'); % Hp
        Obs_interp(3,iperiod, :, iH) = interp2(X, Y, Hr, phi_gm, theta_gm,  'linear'); % Hr
        Obs_interp(4,iperiod, :, iH) = interp2(X, Y, Et, phi_gm, theta_gm,  'linear'); % Et
        Obs_interp(5,iperiod, :, iH) = interp2(X, Y, Ep, phi_gm, theta_gm,  'linear'); % Ep
        
        PeriodArray(iperiod, iH) = Period;
        
        
    end;
    
    iperiod = iperiod + 1;
    
end;


    
% Now write the data to a file


%Electric
save(['/Users/manojnair/projects/GIC/Extracted_GNU/ELECTRIC/ELEC_GNU_Sq' datestr(mean(fday),'mmmYY') '_' Output_directory_name ], ... 
    'Obs_interp', 'NPeriods', 'NHarmonics', 'names', 'phi_gg', 'theta_gg', 'PeriodArray');

%Magnetic

% names = name;
% % 
% save(['/Users/manojnair/projects/GIC/Extracted_GNU/MAGNETIC/MAG_GNU_Sq' datestr(mean(fday),'mmmYY') '_' Output_directory_name ], ... 
%     'Obs_interp', 'NPeriods', 'NHarmonics', 'names', 'phi_gg', 'theta_gg', 'PeriodArray');
% 










